1 /* Copyright 2002-2016 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.semianalytical.dsst.utilities;
18
19 import java.io.Serializable;
20 import java.util.ArrayList;
21
22 import org.hipparchus.analysis.interpolation.HermiteInterpolator;
23 import org.hipparchus.util.FastMath;
24 import org.orekit.time.AbsoluteDate;
25
26 /** Interpolated short periodics coefficients.
27 * <p>
28 * Representation of a coefficient that need to be interpolated over time.
29 * </p><p>
30 * The short periodics coefficients can be interpolated for faster computation.
31 * This class stores computed values of the coefficients through the method
32 * {@link #addGridPoint} and gives an interpolated result through the method
33 * {@link #value}.
34 * </p>
35 * @author Nicolas Bernard
36 *
37 */
38 public class ShortPeriodicsInterpolatedCoefficient implements Serializable {
39
40 /** Serializable UID. */
41 private static final long serialVersionUID = 20160319L;
42
43 /**Values of the already computed coefficients.*/
44 private ArrayList<double[]> values;
45
46 /**Grid points.*/
47 private ArrayList<AbsoluteDate> abscissae;
48
49 /**Number of points used in the interpolation.*/
50 private int interpolationPoints;
51
52 /**Index of the latest closest neighbor.*/
53 private int latestClosestNeighbor;
54
55 /**Simple constructor.
56 * @param interpolationPoints number of points used in the interpolation
57 */
58 public ShortPeriodicsInterpolatedCoefficient(final int interpolationPoints) {
59 this.interpolationPoints = interpolationPoints;
60 this.abscissae = new ArrayList<AbsoluteDate>();
61 this.values = new ArrayList<double[]>();
62 this.latestClosestNeighbor = 0;
63 }
64
65 /**Compute the value of the coefficient.
66 * @param date date at which the coefficient should be computed
67 * @return value of the coefficient
68 */
69 public double[] value(final AbsoluteDate date) {
70 //Get the closest points from the input date
71 final int[] neighbors = getNeighborsIndices(date);
72
73 //Creation and set up of the interpolator
74 final HermiteInterpolator interpolator = new HermiteInterpolator();
75 for (int i : neighbors) {
76 interpolator.addSamplePoint(abscissae.get(i).durationFrom(date), values.get(i));
77 }
78
79 //interpolation
80 return interpolator.value(0.0);
81
82 }
83
84 /**Find the closest available points from the specified date.
85 * @param date date of interest
86 * @return indices corresponding to the closest points on the time scale
87 */
88 private int[] getNeighborsIndices(final AbsoluteDate date) {
89 final int sizeofNeighborhood = FastMath.min(interpolationPoints, abscissae.size());
90 final int[] neighborsIndices = new int[sizeofNeighborhood];
91
92 //If the size of the complete sample is less than
93 //the desired number of interpolation points,
94 //then the entire sample is considered as the neighborhood
95 if (interpolationPoints >= abscissae.size()) {
96 for (int i = 0; i < sizeofNeighborhood; i++) {
97 neighborsIndices[i] = i;
98 }
99 } else {
100 // get indices around closest neighbor
101 int inf = getClosestNeighbor(date);
102 int sup = inf + 1;
103
104 while (sup - inf < interpolationPoints) {
105 if (inf == 0) { //This means that we have reached the earliest date
106 sup++;
107 } else if (sup >= abscissae.size()) { //This means that we have reached the latest date
108 inf--;
109 } else { //the choice is made between the two next neighbors
110 final double lowerNeighborDistance = FastMath.abs(abscissae.get(inf - 1).durationFrom(date));
111 final double upperNeighborDistance = FastMath.abs(abscissae.get(sup).durationFrom(date));
112
113 if (lowerNeighborDistance <= upperNeighborDistance) {
114 inf--;
115 } else {
116 sup++;
117 }
118 }
119 }
120
121 for (int i = 0; i < interpolationPoints; ++i) {
122 neighborsIndices[i] = inf + i;
123 }
124
125 }
126
127 return neighborsIndices;
128 }
129
130 /**Find the closest point from a specific date amongst the available points.
131 * @param date date of interest
132 * @return index of the closest abscissa from the date of interest
133 */
134 private int getClosestNeighbor(final AbsoluteDate date) {
135 //the starting point is the latest result of a call to this method.
136 //Indeed, as this class is meant to be called during an integration process
137 //with an input date evolving often continuously in time, there is a high
138 //probability that the result will be the same as for last call of
139 //this method.
140 int closestNeighbor = latestClosestNeighbor;
141
142 //case where the date is before the available points
143 if (date.compareTo(abscissae.get(0)) <= 0) {
144 closestNeighbor = 0;
145 }
146 //case where the date is after the available points
147 else if (date.compareTo(abscissae.get(abscissae.size() - 1)) >= 0) {
148 closestNeighbor = abscissae.size() - 1;
149 }
150 //general case: one is looking for the two consecutives entries that surround the input date
151 //then one choose the closest one
152 else {
153 int lowerBorder = latestClosestNeighbor;
154 int upperBorder = latestClosestNeighbor;
155
156 final int searchDirection = date.compareTo(abscissae.get(latestClosestNeighbor));
157 if (searchDirection > 0) {
158 upperBorder++;
159 while (date.compareTo(abscissae.get(upperBorder)) > 0) {
160 upperBorder++;
161 lowerBorder++;
162 }
163 }
164 else {
165 lowerBorder--;
166 while (date.compareTo(abscissae.get(lowerBorder)) < 0) {
167 upperBorder--;
168 lowerBorder--;
169 }
170 }
171
172 final double lowerDistance = FastMath.abs(date.durationFrom(abscissae.get(lowerBorder)));
173 final double upperDistance = FastMath.abs(date.durationFrom(abscissae.get(upperBorder)));
174
175 closestNeighbor = (lowerDistance < upperDistance) ? lowerBorder : upperBorder;
176 }
177
178 //The result is stored in order to speed up the next call to the function
179 //Indeed, it is highly likely that the requested result will be the same
180 this.latestClosestNeighbor = closestNeighbor;
181 return closestNeighbor;
182 }
183
184 /** Clear the recorded values from the interpolation grid.
185 */
186 public void clearHistory() {
187 abscissae.clear();
188 values.clear();
189 }
190
191 /** Add a point to the interpolation grid.
192 * @param date abscissa of the point
193 * @param value value of the element
194 */
195 public void addGridPoint(final AbsoluteDate date, final double[] value) {
196 //If the grid is empty, the value is directly added to both arrays
197 if (abscissae.isEmpty()) {
198 abscissae.add(date);
199 values.add(value);
200 }
201 //If the grid already contains this point, only its value is changed
202 else if (abscissae.contains(date)) {
203 values.set(abscissae.indexOf(date), value);
204 }
205 //If the grid does not contain this point, the position of the point
206 //in the grid is computed first
207 else {
208 final int closestNeighbor = getClosestNeighbor(date);
209 final int index = (date.compareTo(abscissae.get(closestNeighbor)) < 0) ? closestNeighbor : closestNeighbor + 1;
210 abscissae.add(index, date);
211 values.add(index, value);
212 }
213 }
214 }